#include <iostream>
#include <fstream>
#include <sstream>
#include <cstring>
#include <cmath>
#include "TFile.h"
#include "TH1D.h"
#include "covariance.h"
#include "matrix_inverse.h"


using namespace std;

//invert matrix and return determinant
double Invert_Matrix(double* matrix, int rank)
{
  MatrixClass<double> M(rank, rank);
  
  //dump matrix entries to external class
  for(int i=0; i<rank; i++)
  for(int j=0; j<rank; j++)
    M.setvalue(i,j,matrix[i*rank+j]);
  
  //invert matrix
  double matrix_determinant=M.invert();
  assert(matrix_determinant!=0);
  
  bool state;
  //dump values back to matrix
  for(int i=0; i<rank; i++)
  for(int j=0; j<rank; j++)
    {
      M.getvalue(i,j,matrix[i*rank+j], state);
      //make sure there is no error
    }

  return matrix_determinant;
}





int main(int argc, char** argv)
{
  if(argc<2+1)
    {
      cout<<"Usage: ./Fisher.exe signal.txt, input1.txt, input2.txt ... "<<endl;
      return 0;
    }

  cout<<"number of input files: "<<argc-2<<endl;
  cout<<"signal: "<<argv[1]<<endl;
  cout<<"first input: "<<argv[2]<<endl;
  
  //grab a list of covariance matrix
  const int rank=4;

  Covariance* sig=new Covariance(rank);
  sig->Read(argv[1]);


  cout<<"Processing file: "<<argv[2]<<endl;
  Covariance* bkg=new Covariance(rank);
  bkg->Read(argv[2]);

  //loop over input
  for(int i=3; i<argc; i++)
    {
      Covariance temp_bkg(rank);
      temp_bkg.Read(argv[i]);

      (*bkg) += temp_bkg;
    }


  

  cout<<"total covariance: "<<endl;
  cout<<*bkg<<endl;

  //now computing the fisher matrix
  double fisher[rank][rank]={};
  double mean[rank]={};

  cout<<"signal mean: "<<endl;
  for(int i=0; i<rank; i++)
    {
      cout<<sig->mean(i);
      if(i<rank-1)
	cout<<", ";
    }
  cout<<endl<<endl;

  cout<<"background mean: "<<endl;
  for(int i=0; i<rank; i++)
    {
      cout<<bkg->mean(i);
      if(i<rank-1)
	cout<<", ";
    }
  cout<<endl;


  cout<<"fisher matrix: "<<endl;
  for(int i=0; i<rank; i++)
    {
      mean[i]= sig->mean(i) - bkg->mean(i);
      for(int j=0; j<rank; j++)
	{
	  fisher[i][j]= (*sig)(i,j) + (*bkg)(i,j);

	  cout<<fisher[i][j];
	  if(j<rank-1)
	    cout<<", ";
	}
      cout<<endl;
    }
  cout<<endl;

  cout<<"determinant: "<<Invert_Matrix(fisher[0], rank)<<endl;
  cout<<endl;

  cout<<"fisher matrix inverse: "<<endl;
  for(int i=0; i<rank; i++)
    {
      mean[i]= sig->mean(i) - bkg->mean(i);
      for(int j=0; j<rank; j++)
	{
	  cout<<fisher[i][j];
	  if(j<rank-1)
	    cout<<", ";
	}
      cout<<endl;
    }
  cout<<endl;



  double result[rank]={};
  for(int i=0; i<rank; i++)
    {
      for(int j=0; j<rank; j++)
	{
	  double old_var=result[i];

	  result[i]+= fisher[i][j]*mean[j];
	  
	  if(isnan(result[i]))
	    {
	      cout<<endl;
	      cout<<"ERROR: nan"<<endl;
	      cout<<"before: "<<old_var<<endl;
	      cout<<"i,j: "<<i<<", "<<j<<endl;
	      cout<<"fisher[i][j]: "<<fisher[i][j]<<endl;
	      cout<<"mean[j]: "<<mean[j]<<endl;
	      cout<<endl;
	    }
	}
    }

  //compute offset and normalization
  double bkg_fisher=0;
  double sig_fisher=0;

  for(int i=0; i<rank; i++)
    {
      bkg_fisher+=result[i]*bkg->mean(i);
      sig_fisher+=result[i]*sig->mean(i);
    }  
  
  double norm=2/(sig_fisher-bkg_fisher);

  
  cout<<"fisher: "<<endl;
  for(int i=0; i<rank; i++)
    {
      result[i]*=norm;
      cout<<result[i];
      if(i<rank-1)
	cout<<", ";
    }

  cout<<endl;
  cout<<"offset: "<<(sig_fisher + bkg_fisher)/(bkg_fisher - sig_fisher);
  cout<<endl;
  cout<<"relative importance: "<<endl;
  for(int i=0; i<rank; i++)
    {
      cout<<result[i]*(sig->mean(i)-bkg->mean(i));				
	//*((*sig)(i,i)+(*bkg)(i,i));
      if(i<rank-1)
	cout<<", ";
    }

  cout<<endl;

  
    
  
  

  return 1;
}
